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Ballistic annihilation with continuous initial velocity distributions is investigated in the framework 
of Boltzmann equation. The particle density and the rms velocity decay as c ~ t~ a and (v) ~ t~^ , 
with the exponents depending on the initial velocity distribution and the spatial dimension. For 
instance, in one dimension for the uniform initial velocity distribution we find /3 = 0.230472 .... We 
also solve the Boltzmann equation for Maxwell particles and very hard particles in arbitrary spatial 
dimension. These solvable cases provide bounds for the decay exponents of the hard sphere gas. 

PACS numbers: 05.20.Dd, 03.20.+i, 82.20Mj 



Ballistic annihilation is the kinetic process which in- 
volves particles undergoing ballistic motion and anni- 
hilating upon colliding. Ballistic annihilation underlies 
numerous apparently unrelated phenomena, e.g., growth 
and coarsening processes and traffic flows flS-El . In coars- 
ening processes, for instance, domain walls have a nat- 
ural particle interpretation. Despite all that, little is 
known on irreversible processes where the reactants move 
ballistically while the contrasting situation of diffusion- 
controlled processes is well understood || . 

There exist few theoretical results on ballistic annihila- 
tion in one dimension with discrete velocity distributions. 
For the simplest binary velocity distribution, the ballis- 
tic annihilation process has been solved in one dimension 
by Elskens and Frisch ||, see also Refs. some 
analytical results are also available for the ternary veloc- 
ity distribution 0,0] • No solutions have been found for 
continuous initial velocity distributions, although the de- 
cay exponents have been determined numerically Jl2| ? |l3[ | . 
This lack of analytical results is especially striking given 
that ballistic aggregation processes admit exact solutions 
in one dimension for arbitrary initial velocity distribu- 
tions ||,|| . 

In this work, we consider ballistic annihilation with 
continuous isotropic initial velocity distributions in arbi- 
trary dimension. Our analysis is performed in the frame- 
work of Boltzmann equation approach. This scheme in- 
volves an uncontrolled approximation and generally leads 
to erroneous results for ballistic annihilation with discrete 
velocity distributions. For continuous velocity distribu- 
tions, however, the decay exponents found from numer- 
ical integration of the Boltzmann equation E] a re an 
excellent agreement with simulation results |12fll3|| . Fur- 
thermore, for solvable one dimension ballistic aggregation 
processes with continuous velocity distributions 0,01, 
the decay exponents computed from the Boltzmann equa- 
tion are exact; for higher-dimensional ballistic aggrega- 
tion process [0] , the decay exponents determined within 
the Boltzmann framework [0,0 appear to be exact as 
well. Therefore, we believe that the Boltzmann equation 
approach provides exact decay exponents for ballistic an- 
nihilation with continuous velocity distributions. 

For clarity, we start with one-dimensional ballistic an- 



nihilation process. An appropriate Boltzmann equation 
reads [HH 



dP(v,t) 
dt 



P(v,t) / dv'\v-v'\P(v',t). 



(1) 



In the long time limit, the velocity distribution ap- 
proaches a scaling form 



P( V} t) = ^~ Q F(x), with x = vt 13 . 



(2) 



The two basic exponents determine the behavior of the 
particle concentration c(t) and the rms velocity (v): 



Formally, c and (v) are defined as follows 



(3) 



(t) = / dvP(v,t) = t~ a 



dx F(x) 



(v) 2 = c- 1 



dvv 2 P(v,t) = t 



- f-2/9 



f^ > oo dxx 2 F(x) 

rLdxF( X ) ■ 



By inserting the scaling form (|J) into Eq. (|l]) one finds 
a + [3 = 1 and an equation for the scaling function 



2/3 - 1 + (3x 



F'(x) 
F(x) 



dx'\x-x'\F{x'). (4) 



In the following, we always consider isotropic initial ve- 
locity distributions. In one dimension, this requirement 
reads Po(v) — Pq(—v) and it implies the symmetry for 
later times P(v,t) — P(—v,t) and the symmetry of the 
scaling function F(x) = F(—x). 

The large x behavior of the scaled velocity distribution 
is found by noting that in this region the integral on the 
right-hand side of Eq. (|j) simplifies to Cx, where C is 
the normalization constant, C = dx F(x). Solving 
the resulting differential equation gives 



F(x) ~ e- Cx ' when x - oo. 



(5) 



For a given (3, a solution of Eq. (Q) would not agree 
with the boundary condition at x = which is implied 



1 



by the initial velocity distribution. Thus, we arrive at 
an eigenvalue problem. For instance, we should require 
F'(0) = if the initial velocity distribution is flat near the 
origin. Solving this eigenvalue problem numerically gives 
/? flat = 0.230472 . . ., to be compared with /3 flat » 0.22 
pj and /3 flat « 0.19 |p| found from Monte Carlo sim- 
ulations. Similarly, for the initial velocity distribution 
satisfying Pq(v) ~ \v\ we find Ainear = 0.166649 . . .. 

We now outline an approximate analytical computa- 
tion of the exponent [3. First we note that Eq. (|J) can 
be reduced to an ordinary differential equation after a 
double differentiation. This equation can be further sim- 
plified by the rescaling x — > Xy/]3 that eliminates the j3 
factor. The governing equation then reads 



xf"' + 2f" = 2exp(-/), 



(6) 



with f(x) = — \nF(x). We shall also use the relation 

dxxF(x) = (2/3)- 1 - 1, (7) 



which will play the role of a normalization condition. 
Equation (|7) is just Eq. ([|) at x = 0. 

Although Eq. (^|) cannot be solved exactly, an approxi- 
mate solution fx can be found by replacing the right-hand 
side by e~^°, where /o is a reasonable approximation for 
/. In principle, this approximation scheme can be re- 
peated again starting from /i, and should finally lead to 
the exact form for /. As the starting point, we choose 



\ v 2v A J 



-x I V 



(8) 



The function F (x) is constructed in such a way, that the 
small x expansion is correct up to the second order. Ad- 
ditionally, F (x) exhibits an exponential decay for large 
x, in agreement with the exact asymptotic behavior (|^). 
The parameter v is yet to be determined. Replacing now 
the right-hand side of Eq. (||) by Fq (x) from Eq. (g) and 
solving the resulting linear differential equation gives 

h{x) = 2^(3 -v 2 )x + !/ 2 (7-5z/ 2 ) {l-e-'f v \ 
- v{l-v 2 )xe- x ' v 



6^(2-^) l" V lz^k± 



dt. 



(9) 



The constant v can now be calculated self-consistently 
by imposing the constraint 



/•OO 

= 2 1 F(x) dx, 
Jo 



(10) 



which can be obtained by integrating Eq. (||). Plugging 
F appr = e ~ h ^ into Eq. (0) gives v~^ T = 2.67156 to be 
compared to the exact numerical value v~ x = 2.65826 . . .. 
Finally, Eq. @ leads to /3 apP r = 0.22898, in good agree- 
ment with the exact value obtained above. In Fig. 1, we 



plot the exact numerical scaling function F and the ap- 
proximate solution -Fappr = e^ 1 ^) with f\{x) given by 
Eq. (^). The approximate solution is extremely good as 
the relative error is always less than 0.5%. 




FIG. 1. Plot on a semi- log scale of the scaling function F 
and its approximation given by Eq. (|^). The insert shows 
both functions in normal scale. Both scaling functions are 
indistinguishable as the relative error is less than 0.5% 

In higher dimensions, the Boltzmann equation reads 
<9P(v, t) 



dt 



-P(v,t) J dw|v-w|P(w,f). 



(11) 



For isotropic initial velocity distributions, the appropri- 
ate scaling variable is x — vt" with v — |v|, and the 
scaling form is P(y,t) = t °~ a F(x). Plugging this scal- 
ing form into Eq. ( |ll|) and using the spherical coordinates 
to simplify the collision integral we obtain 



(d+l)/3 -l+/3x 



F'(x) 
F(x) 



(12) 



= -n d _ 1 dyy d - 1 E d {x,y)F{y). 
Jo 

Here Qd-i is the surface area of the unit sphere in d— 1 
dimensions and 



Ed{x,y) 



> (sin ( 



\d-2 



\/ x 1 + y 2 — 2xy cos 4>- (13) 



In two dimensions, Ei(x^ y) = 2(x + y)E(k) 1 where E{k) 
is the complete elliptic integral of the second kind with 



modulus k 



x+y 



Similarly, in an arbitrary even di- 



mension Ed(x,y) can be expressed via elliptic integrals. 
In odd dimensions, Ed{x,y) can be expressed in terms 
of elementary functions. In the most interesting three 
dimensional case, one reduces Eq. (12) to 



4/3 - 1 + j3x 



2vr 



F'(x) 
F(x) 

i 

dyy 



(14) 



(x + yf 



F(y)- 
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One should solve Eq. ( |12| ) subject to an appropriate 
boundary condition at the origin. Overall, the task re- 
duces again to the eigenvalue problem. Note that the ap- 
proximation scheme explicitly presented above in d = 1 
could be equally applied in higher dimensions. 

Thus, the decay exponents can be determined with 
arbitrarily high precision in arbitrary spatial dimension 
for arbitrary isotropic initial velocity distribution. This 
method, however, does not solve the Boltzmann equa- 
tion (|l|). We now provide two approximate solutions to 
the Boltzmann equation for arbitrary initial conditions 
and arbitrary spatial dimension d. These solutions are 
in fact exact solutions of Boltzmann equations with col- 
lision kernels a(g) different from a(g) = g characterizing 
the hard sphere gas (here g = |v — w| is the relative ve- 
locity). Note that such collision kernels naturally arise in 
kinetic theory of interacting particles |l£| l . In the present 
context, the generalized Boltzmann equation reads 



dP(v,t) 
dt 



P(v,t) / dwa(g)P(w,i). 



(15) 



Let us compare dimensions of the left and right-hand 
side of Eq. (p"5|). The velocity distribution has dimen- 
sion [P] = T/L d+1 which implies [a] = L d /T; therefore 
[o~(g)/g] = U 1 ^ 1 . The remaining quantity with dimen- 
sion of length, the "interaction" radius, should be ex- 
tracted from the collision process. For hard sphere gas, 
the relevant interaction radius is simply the geometri- 
cal radius a of the spheres, so <r(g) — ga^ 1 . The con- 
stant factor a d_1 can be absorbed into the time vari- 
able - this is what we have done in Eq. (|TT|) . For parti- 
cles interacting through a two-body power law potential, 
U(r) oc r~ n , the energy conservation implies g 2 ~ r~ n . 

Thus, <r(g) ~ gr^ 1 ~ g v with v = 1 - ^T^- Thc hard 
sphere gas (y = 1) is recovered for n = oo. The veloc- 
ity independent kernel {y = 0), the so-called "Maxwell" 
gas, arises when particles interact through the power 
law potential with the exponent n = 2{d — 1). When 
v > 1, the interaction is "harder" than in the hard sphere 
gas (though such behavior does not arise from a sim- 
ple power-law interaction potential). One particularly 
tractable model corresponds to a(g) = g 2 , the so-called 
gas of very hard particles . 

To provide a faithful analog of the original hard sphere 
gas, we replace a(g) — g by (v) in the Maxwell case and 
by the factor g 2 /(v) for the very hard particles. Hence for 
the Maxwell gas the Boltzmann equation (|l|) becomes 



dP(v, t) 
dt 



= -P(v,t)c(t)(v(t)). 



(16) 



Thus we get effectively non-interacting particles as dif- 
ferent velocities remain uncoupled. Solving (|l6|) yields 



P(v,t)=c(t)P (v), c(t) = 



1 



l + (v)t 



(17) 



The moments of the velocity do not change with time, 
e.g., (v(t)) = (v)o, and thence a = 1 and = 0. 



More interesting results are found for the very hard 
particles. The corresponding Boltzmann equation, i.e., 
Eq. ( |l5| ) with a(g) = g 2 /(v), can be simplified by ab- 
sorbing the (v}~ 1 factor into the time variable, 



dt' 



o (v(t')) , 



(18) 



and reducing the collision integral into a combination of 
the moments Mj(r) = J dw P(w,t) of the velocity 
distribution. Thc Boltzmann equation becomes 



dPjyyr) 
dr 



-P(v,t)[v 2 M (t)+M 2 (t)], (19) 



and easily solved to give 

P(v,t) = P (v)e- 1,2io(r) - i2(r) . 



(20) 



Here L 3 {t) = f Q dr' Mj(t'). _ 

To derive explicit results, it is natural to consider ini- 
tial velocity distributions algebraic near the origin. To 
simplify algebra, we specifically choose 



Po(v) = 



2v»e~ v 



a d T [(/x + d)/2] 



(21) 



where n d = 2n d / 2 /T(d/2) is the surface area of thc unit 
sphere in d dimensions, and the prefactor in ( ^l| ) is cho- 
sen to set the initial density to unity. Combining J2(i| ) 
and (Ell) we can explicitly compute Mq and M^- Using 



dL 



then Mj — closes the problem. We find 



dL 
dr 



dL? 



Solving these equations subject to Lq(0) = L 2 (0) = 
yields 

L = [1 + (fi + d+ l) T ]S+3+T - l, 

L -27^Tij ln[1 + ^ + rf + 1) ^ 

Now we compute c = Mq = ^ tt to find the density, 

c = [1 + (fj, + d + l) T }~~rf^ , (22) 
and (v) 2 = ^ = ^(l + L )-\ 



[1 + (fi + d + 1)t] 2 (tl+ d+i) 



(23) 



By inserting Eq. ( p3| ) into Eq. (|T^) we can express r via 
the original time variable t. Finally, we arrive at 



2(p + d) 
2{n + d) + l' 



2(/i + d) + l 



(24) 
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Note that the exponent relation a + (3 = 1 which is valid 
for any fi and d. This sum rule immediately follows from 
an elementary mean-free path argument: a d ~ 1 c(v)t ~ 1. 

The above exact values of the exponents in the two 
solvable limits appear to provide the strict bounds for 
the hard sphere case: 



2(p + d) 
20 + d) + 1 



< a < 1, < (3 < 



1 



2(/i + d) + 1 ' 



(25) 



+ |v'| K ) (the most natural case 
For k > 0, the decay exponents 

2n 



These bounds are fair for small d and they get more and 
more stringent as the spatial dimension increases. 

The above method of solving the Boltzmann equa- 
tion can be adapted to the more general collision kernels 

cr(v,v') = (w) 1_K (|v| K 

corresponds to k = 1). 

are given by Eqs. d2~|). The collision kernels a(g) = g 
constitute a more perspective generalization. One could 
try to solve the Boltzmann equation when n is integer 
and then perform an analytic continuation to n = 1/2 
corresponding to the hard sphere gas. Besides the cases 
of Maxwell and very hard particles (n = and 1, re- 
spectively), it is possible to work out the case of n = 2. 
Unfortunately, we have not succeeded beyond that. 

In the limit d — > oo, the Boltzmann equation ( jll| ) be- 
comes tractable. First of all, the collision kernel simplifies 
to \J v 2 + w 2 as (different) vectors are orthogonal in infi- 
nite dimensions. The scaled Boltzmann equation reads 



F' 

l-B-px— = fl d 



dyy d -Wx 2 + y 2 F(y), (26) 



with B = (d+ 1)(3. The bounds of Eq. (||) lead to j3 -> 
but remain non-trivial for B, < B < 1/2. The right- 
hand side of Eq. ( p6| ) is computed by the saddle point 
technique to find C y 1 x 1 + y 2 , where C is the normaliza- 
tion factor C — fid Jq 00 dyy^ 1 F(y), and y* is the saddle 
point which is found from y* F'{y*)/F{y*) + d — 1 = 0. 
Hence, at the saddle point Eq. ( p6| ) gives 1 = Cy«y/2. 
Near the origin, F(x) ~ x' 1 and thus 0xF'/F = — > 0, 
so we find 1 — B = Cy* = 1/ y/2. Thus we arrive at the 
universal asymptotics, 



i - 4= ) '/ 

V2y 



when <i 



oo, 



(27) 



independently on fi. 

In this work, we reduced the determination of the de- 
cay exponents for ballistic annihilation to an eigenvalue 
problem. We found that the exponents have very non- 
trivial values even for the simplest initial velocity distri- 
butions. Our approach manifestly demonstrates that the 
decay exponents are affected only by the spatial dimen- 
sion d and by the exponent fi characterizing the initial 



velocity distribution in the |v| — > limit: Pq(v) ~ |v| M . 
We also solved the Boltzmann equation for the Maxwell 
particles and very hard particles in arbitrary spatial di- 
mension. For the hard sphere gas, we found the asymp- 
totic behavior of the exponents in the d — > oo limit. 

Interesting discussions with Eli Ben-Nairn and Satya Ma- 
jumdar are gratefully acknowledged. PLK also acknowl- 
edges support from NSF, ARO, and CNRS. 
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